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ABSTRACT 

Starting with a mathematical boundary value problem for the magnetic vector potential in an axisym- 
mctric cylindrical coordinate system, we derive a general solution for any arbitrary current distribution 
using the method of Green's functions. We use this to derive an analytic form for an hourglass mag- 
netic field pattern created by electrical currents that are concentrated near (but not confined within) 
the equatorial plane of a cylindrical coordinate system. Our solution is not characterized by a cusp 
at the equatorial plane, as in previous solutions based on a current sheet. The pattern we derive pro- 
vides a very good fit to hourglass magnetic field patterns emerging from three-dimensional numerical 
simulations of core formation, and can in principle be used for source-fitting of observed magnetic 
hourglass patterns. 

Subject headings: stars: formation — ISM: magnetic fields — ISM: individual objects (NGC 1333 
IRAS 4A, G31.41-I-0.31) — methods: analytical 



1. INTRODUCTION 

An hourglass-shaped magnetic field is a natural out- 
come of gravitationally-driven contraction of a plasma 
that is embedded in a large-scale magnetic field. Self- 
inductance and high conductivity of the plasma leads to a 
flux-freezing effect that drags the magnetic field lines in- 
ward in the dense contracted regions. Hourglass patterns 
are evident in theoretical calculations of magnetohydro- 
static equilibria (Mouschovias 1970) and in dynamical 
calculations of collapsing magnetized molecular cloud 
cores (e.g., Kudoh & Basu 2008, 2011). The first un- 
ambiguous observed hourglass pattern in a star-forming 
region was in NGC 1333 IRAS 4A (Girart ct al. 2006). 
This was measured by submillimeter dust emission that 
is polarized perpendicular to the direction of the ambient 
magnetic field. Subsequent dust polarization measure- 
ments also show an hourglass pattern in the high-mass 
star-forming region G31.41-h0.31 (Girart et al. 2009). 

These measurements confirm the general scenario of 
gravitational collapse in regions with large-scale mag- 
netic fields, and that the turbulent energy in the plasma 
does not dominate the magnetic energy. The ambient 
strength of the magnetic field is however not measur- 
able through the polarization signal, due to uncertain- 
ties in grain properties that lead to alignment and emis- 
sion strength. The shape of the hourglass pattern does 
however hold the potential to constrain the magnetic 
field strength (see Basu et al. (2009) for a discussion of 
this point based on numerical models). The degree of 
pinching of the hourglass contains information about the 
ratio of central to background magnetic field strength. 
If this is combined with information about the ratio of 
densities in these regions, one can estimate the relative 
degree of flux-freezing in the contraction of the star- 
forming region. By this indirect means, one can deter- 
mine the degree of flux-freezing and hence the degree of 
ambipolar diffusion (neutral-ion drift) during the con- 
traction. There arc essentially two limits of core for- 
mation in a magnetized environment: (1) contraction 
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due to neutral-ion drift (with low amounts of inward 
field-line dragging) in a magnetically-dominated envi- 
ronment, and (2) gravitationally-dominated contraction 
with significant inward field-line dragging. The former 
occurs on the ambipolar diffusion time scale and the lat- 
ter on the (shorter) dynamical time scale. These are 
well illustrated through linear analyses of gravitational 
instability in magnetized media (Ciolek & Basu 2006; 
Mouschovias et al. 2011; Bailey & Basu 2012). There is 
a continuum of possibilities between the above two lim- 
its, with increasing amounts of pinching of the hourglass 
while progressing from the first to second limit. 

In this paper, we find an analytic form for an hour- 
glass magnetic field pattern, extending previous studies 
(Mestcl & Ray 1985; Barker & Mestel 1990) that have 
calculated patterns based on a current sheet in the equa- 
torial plane of a cylindrical coordinate system. Unlike 
those solutions, the pattern derived here does not have a 
cusp at the midplane, and is more generally able to fit re- 
sults of three-dimensional numerical simulations. The in- 
tent of this work is to provide analytic formulae that can 
be used relatively easily for observational source-fitting 
and characterizing properties of the magnetic field, e.g., 
the ratio of central to background field strength. Our 
model is formally suited to fit the hourglass magnetic 
field pattern of a prestellar core, however it could also be 
applied to protostellar cores in some cases, as discussed 
in Section 4. 

A theoretical hourglass pattern can also be used to 
generate synthetic polarization emission maps, as done 
recently by Kataoka et al. (2012) using numerical models 
of hourglass fields. Their numerical models also include a 
toroidal magnetic field that is generated by rotation, and 
they point out the potential of source-fitting in order to 
gain information about magnetic field strength, rotation 
rate, and evolutionary state. 

2. MODEL 

Following Mestel & Ray (1985), we assume that the 
magnetic field threading a flattened star-forming core 
(hereafter referred to as "the disk") can be considered 
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as a local distortion of an initially uniform background 
field Bq. We introduce cylindrical coordinates {r,ip,z), 
with the origin at the center of the protostellar disk and 
the z-axis parallel to the background field. Thus, the 
background field is 



(1) 



where Bq is a constant. We express the net magnetic 
field B as 

B = Bd + B„, (2) 

where Bd represents the local magnetic field generated 
by the disk currents. 

We require that the magnetic field through the disk is 
axisymmetric and thus define 



d_ 

dip 



0. 



(3) 



In order to produce a simple hourglass morphology with- 
out any azimuthal twisting in the field lines, we also re- 
quire that B is purely poloidal, that is 

B = Brr + B,z. (4) 

Since B, and thus B^, must be solenoidal, we write 

Bd = VxA, (5) 

where A is a purely toroidal Coulomb gauge vector po- 
tential, given by 



A = Aip = — (p . 

r 



(6) 



The function $(r, z) is the poloidal magnetic flux (di- 
vided by a factor of 27r) due to Bd, and can be expressed 



as 



$ = / Bd- z r'dr' . 
Jo 



(7) 



Hereafter we will refer to $ simply as "the flux." For a 
derivation of Eq. (6), see Kulsrud (2004). 

In general, using the pre-Maxwell form of Ampere's law 
(expressed in Gaussian cgs units), the Coulomb gauge 
vector potential A produced by a volume current density 
j satisfies 



^2 A 47r . 

c 

V- A = 0. 



(8) 
(9) 



Given our expression for A in Eq. (6), it follows that 
the volume current density in the disk must be purely 
toroidal as well, and we write 

j = j<P- (10) 
Thus, Eq. (8) simplifies to 

A A-K 

We treat r and z as ordinary Cartesian coordinates, 
and introduce the operator V with respect to these vari- 
ables: 



(11) 
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Multiplying Eq. (11) by a factor of r, we express the 
PDE in the form 

LA = -V-(rVAj+^ = ^rj, (13) 

where L is a self-adjoint elliptic operator. This greatly 
simplifies the task of finding the Green's function for this 
problem. 

We impose the boundary conditions that B must relax 
to the background field Bq as z — >■ ±oo, and likewise for 
r > R, where R is the disk radius. Thus, we require that 
the locally generated flux satisfies 



*(0,z) = 0, 
<i>(i?,z)=0, 
lim $(r, z) = 0. 

^— >-±oo 



(14) 
(15) 
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Using Eq. (6), we translate the preceding conditions on 
$ into equivalent conditions on A: 



A{0, z) bounded, 
A{R,z) = 0, 
lim A(r, z) = . 

2— >-±00 



(17) 
(18) 
(19) 



Given these conditions, we solve Eq. (13) for A on the 
domain D = { (r, z) e | r e (0, i?), zG (-00,00)}. 
We do so by finding the Green's function G(^, 77, r, z) for 
the problem. Thus, we first solve the associated PDE 



LG{tii,T,z)^6{r-05{z~il) 



(20) 



for a point source located at (r, z) = (^, rj) (z D with 
identical boundary conditions (17) through (19) imposed 
on G{^,r],r, z). Since the coefficients in the operator L 
and the boundary conditions have no explicit dependence 
on z, we write Eq. (20) in the form 



LG = - 



+ LrG = S{r-OS{z-v), 



where the self-adjoint operator Lr is defined by 

L G^-— (r—\ +- 
^ dr \ dr I r 



(21) 



(22) 



We express the Green's function in terms of the eigen- 
function expansion 



G{^,'q,r,z) 



m—l 



Cm(^,f?, 2)?/'„i(r) . 



(23) 



The eigenfunctions ipmif) corresponding to eigenval- 
ues Am are obtained from the following Sturm-Liouville 
eigenvalue problem associated with the operator Lr'. 

Lrlpm = XnTll^m , (24) 

for r G (0, R) with the homogeneous boundary conditions 

%pm{0) bounded, (25) 
V'm(i?) = 0. (26) 

This problem amounts to solving Bessel's differential 
equation with the prescribed boundary conditions. The 
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normalized cigcnfunctions are thus well known (Zauderer 
2006), and are given by 



V2 JiiVK^r) 

R \j2{Vk:r)\ ' 



(27) 



where we use the notation Jq to denote Bessel functions 
of the first kind of order a. The corresponding eigenval- 
ues are 
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where am,i is the m*'' positive root of Ji{x) and a„i^i < 
ttm+i,!- The eigenfunctions are orthonormal with respect 
to the r- weighted inner product. Explicitly, 



(V'm, V'ri 



ipm{r)tpn{r) r dr = 5„ 



(29) 



where Smn denotes the Kronecker delta. 

Having expressed the Green's function in terms of the 
eigenfunction expansion in Eq. (23), we take advantage 
of the orthonormality condition by multiplying Eq. (21) 
by '0T1 {i") and integrating with respect to r over the inter- 
val (0,i?). This yields an ordinary differential equation 
for the expansion coefficients Cni^,!], z): 



+ XnCn{^,V,z)^S{z-r^)MO- (30) 



This in turn implies that the coefficients can be expressed 
in the separable form 



C'ni^,T],z) = K„(77,z)?/'«(6) : 



(31) 



where Kn{r], z) is a Green's function determined from the 
problem 



d'^Kn{ri,z) 
dz^ 



+ A„A'„(7/,z) = <5(z- 77). (32) 



Using the standard Fourier transform method (Zauderer 
200G), we find the solution to be 



2\/A^ 



(33) 



We now combine Eq. (23), (27), (31), and (33) to get 
an explicit expression for the Green's function: 



G(e,r,,r,z)= ^ 



MVK^r) Ji{VK:0 e-v^l--"l 



m — 1 



(34) 

In terms of the Green's function, the general solution to 
the PDE (13) with the associated boundary conditions 
(17) through (19) on the domain D is 
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G{^,r,,r,z)j{^,r^)^d^df^. (35) 



-oo 



For a derivation of the general form of Eq. (35), see 
Zauderer (2000). 

We assume that the current density has the separable 
form 

jir,z) = f{r)giz). (36) 



In this paper, we assume that the ^-dependence of the 
current can be modelled by the Gaussian function 



(37) 



where /i is a free parameter. This choice results in 
a model with smooth field lines. In contrast, setting 
g{z) = 8{z) would yield a thin disk model related to those 
presented by Mestcl Ray (1985) and Barker fc Mcstel 
(1990), which exhibit cusps in the field lines at the mid- 
plane. 

We insert our expressions in Eq. (34) and (36) into Eq. 
(35), which we can then express as 



where we have defined the coefficients 



(38) 



cR^V>^[J2{VK^R)\' 



(39) 



The extra exponential and constants which appear in Eq. 
(38) cancel with their counterparts in the coefficients k^- 
We introduced these factors in order to simplify the form 
of the final solution. 

The integral with respect to rj in Eq. (38) is generally 
much easier to compute by expressing it as 



/(t?) e-^l^-^ldT^ == C {g{z - u) + g{z + u)} 



(40) 

where £ denotes the unilateral Lapl ace transform oper- 
ator from the u domain to the ^/Am domain, defined by 



oo 

£61 = y e{u)e-''^du. 



(41) 



Inserting our choice of g{z) into Eq. (38) and integrating 
leads to the final solution 



A{r,z) = y^fc,nJi(vX^r) 



m—l 
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A,„ h z 



(42) 



The magnetic field components, including the local and 
background contribution, are given by the expressions 
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Figure 1. Plot of the model magnetic field lines. The model 
parameters are fci = 0.9549, fc2 = 0.4608, fca = 0.6320, R = 3.8918, 
h = 0.3257, and Bq = 3.3118, taken from the best fit presented 
in Fig. 2. The values of R and h are measured in the length unit 
adopted by Kudoli Basu (2011). The line spacing is inversely 
proportional to the field strength. 

which yield the cxphcit solutions 



erfc 



Br — ^ ^ fcm (•y/Am r) 



erfc 




(45) 



711 — 1 



erfc 



Am r) 




(46) 



These expressions are valid for (r, z) £ D. For r > R, 
the magnetic field is equal to the background value Bq. 

Rather than assuming a particular model for f(r) and 
computing the integral in Eq. (39) directly, we treat 
the coefficients km as unknown fitting parameters. This 
approach requires that the series is truncated. Fitting 
the truncated series corresponds to fitting longer wave- 
length variations in the data, since the eigenvalues Am 
arc related to the inverse wavelength of the correspond- 
ing eigenfunctions tpm ('') i and smaller eigenvalues are in- 
dexed by smaller m- values. 



3. APPLICATION 



3.1. 



Sample fit to simulation 

The fitting parameters involved in our model are R, h, 
Bq and the coefficients fcm- We found that taking only 
the first three terms in the series for Br and Bz gave 
good results. 

To test our model, we fitted magnetic field data 
from the simulation labeled VO that is presented by 
Kudoli fc Basu (2011). The data was produced by solv- 
ing the three-dimensional MHD equations, including am- 




Figure 2. Plot of the magnetic field lines for the VO simulation 
(red lines) and the best fit model (blue lines) in the plane y = 0. 
The optimized fitting parameters are fci = 0.9549, = 0.4608, 
fcg = 0.6320, R = 3.8918, h = 0.3257, and Bq = 3.3118. The 
spacing of the lines is determined by the simulation grid and not 
the field strength. 

bipolar diffusion and self-gravity. For a complete discus- 
sion, see Kudoli fc Basu (2008). The hourglass magnetic 
field in the VO simulation exhibits remarkable cylindri- 
cal symmetry about the z-axis, and has a negligible az- 
imuthal component, making it an ideal data set for our 
model. 

We took a slice of the data in ?/ = plane and fitted our 
model in the least squares sense. That is, we minimized 
the sum 



(47) 



where Bvo{xi , Zi) denotes the x- and z-components of 
the simulation magnetic field at the point {xi , z^) in the 
plane y = 0, and B(xi , Zi) denotes the model magnetic 
field. We used MATLAB's sequential quadratic program- 
ming algorithm to perform the optimization. 

The field lines generated using the resulting fitting pa- 
rameters are plotted in Fig. 1. The model shows very 
good agreement with the simulation field line morphol- 
ogy, as illustrated in Fig. 2. 

3.2. Fitting polarization data 

Fitting polarization data requires a different ap- 
proach. Polarization maps, such as the ones published 
by Girart et al. (200G, 2009), only convey the morphol- 
ogy of the field lines, but not the field strength. Thus, 
fitting can be performed by minimizing the angle be- 
tween the model magnetic field vectors and the vectors 
obtained from the polarization data. 

The shape of the field lines in our model depends only 
on the ratio B /Bq, regardless of the specific value of the 
background field strength Bq. Thus, when fitting po- 
larization data, Bq should be measured independently 
rather than being determined as a fitting parameter. We 
should instead fit B/Bq to obtain the parameters i?, h, 
and km/Bo. Then, provided an estimate for the back- 
ground field strength, we can calculate the actual B- 
values in the region of interest. 
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Figure 3. Plot of the 2-componcnt of the model magnetic field 
at 2 = 0. The fitting parameters are k\ = 0.9549, k2 = 0.4608, 
fcs = 0.6320, R = 3.8918, h = 0.3257, and Bq = 3.3118. 
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Figure 4. Plot of the r-component of the model magnetic field at 
z = 0.26 iJ. The fitting parameters are fci = 0.9549, fc2 = 0.4608, 
fcs = 0.6320, R = 3.8918, h = 0.3257, and Bo = 3.3118. 

Our fitting function can be compared directly to a 
plane-of-sky polarization map in the manner described 
above. However, a more sophisticated approach would be 
to use a model for thermal dust emission and polarization 
(e.g., Padovani et al. 2012) to combine our axisymmctric 
three-dimensional fitting function with a model of dust 
properties in the region. This would lead to a synthetic 
plane-of-sky polarization map that can be compared with 
the observed polarization map. 

4. CONCLUSIONS AND DISCUSSION 

Our model can be a useful tool for analyzing observa- 
tional data. Provided that an estimate of the background 
field is available, the model can be used to predict the 
central field strength in a cloud. The ratio of central 
to background magnetic field strengths can also be com- 
pared with the ratio of central to background column 
densities in order to determine the degree of flux freez- 
ing- ^ 

Using the Green's function in Eq. (34), it is easy to 
modify our model by using different assumptions about 



the current in the disk. For instance, it would certainly 
be possible to reduce the number of fitting parameters by 
picking a particular model for the r-dependent portion of 
the current, /(r), and evaluating the integral in Eq. (39) 
explicitly. Since the disk has a finite size, /(r) should 
be chosen so that the function goes smoothly to zero at 
r = R. 

Although our model certainly provides an excellent fit 
to the simulation data of Kudoh & Basu (2011), there is 
one drawback that should be discussed. As shown in Fig. 
3, when the truncated series for is fitted to a data set 
the resulting function generally has a small discontinuity 
across r = R, the size of which decreases as z — >■ ±oo. 
This discontinuity occurs because is given by Eq. (43) 
and we assume that the disk has a finite size R. It is im- 
portant to note that the truncated series for S,. is always 
continuous across r = i?, as illustrated in Fig. 4. Con- 
sequently, the discontinuity is in no way apparent when 
plotting the field lines since they are perfectly smooth 
for all r- values. For example, see Fig. 1. 

It is possible to minimize the size of the discontinu- 
ity by introducing the optimization constraint that Bz 
is continuous at (r, z) = (i?, 0). However, we found that 
this additional constraint reduced the overall quality of 
the fit when using three terms in the series. Taking sev- 
eral more terms in the series would likely give better 
results. 

For source-fitting purposes, our Gaussian-disk model 
has two important advantages over the related field line 
model of Barker k Mostd (1990), even though both 
models would yield similar field line patterns at large dis- 
tances from the midplane. That model is infinitesimally 
thin in the vertical direction, leading to a cusp at the 
midplane of the hourglass. Any real star-forming core 
will have a finite vertical thickness, and our fitting pa- 
rameter h can be identified with this physical thickness. 
The Barker & Mcstel (1990) model also has an infinite 
radial extent. Had we also assumed an infinite radial 
extent, then the series over discrete wavelengths in Eq. 
(42) would be replaced by an integral over a continu- 
ous spectrum. The infinite thin-disk model published by 
Barker McstcI (1990), for example, requires numerical 
integration and iteration in order to compute the mag- 
netic field. Our truncated series for the magnetic field, 
on the other hand, is much simpler to evaluate and fit 
to data using basic nonlinear least-squares techniques. 
The fitting parameter R can also be identified with the 
physical radial extent of the star-forming core. 

Our model is formally designed to fit the hourglass 
pattern of a prestellar core, since the profile of B^ is fiat 
(not cusped) at the center, Br goes to zero at r = 0, 
and there is no toroidal field B^. Once a central proto- 
star has formed, the protostellar core will have cusped 
power-laws in B^ and Br at small radii, as well as a 
significant ratio of B^ to B^ (see Machida et al. 2007; 
Dapp ct al. 2012). All of this occurs within the radius of 
the expansion wave that emanates from the center after 
the protostar forms. Our model could be applied to a 
protostellar core if the region of the expansion wave is 
small and/or unresolved. Furthermore, depolarization of 
the inner high density regions may occur due to a de- 
crease of grain alignment efficiency (see Padovani et al. 
2012), so that effects of protostar formation may not be 
easily measurable by polarization maps. 
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